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1 . INTRODUCTION 


This report is concerned v/ith the description of a computer 
program for solving a certain class of optimal control problems 
known as the linear-quadratic problem (LQP) . ,:>I The LQP is 
sometimes referred to as the neighboring optimal guidance 
problem"' ^ even though it is applicable to both guidance and 
control problems. The distinguishing feature of the LQP is that 
it can be solved without iteration (whereas general optimization 
problems usually require iterative numerical techniques). Thus, 

it is useful in the initial portion of guidance and control 
design for determining approximate feedback gains and giving 
insight into the systems. 

The computer program described in Appendix A solves the 
following problem: 

Minimize: J=ixjs f x f +i J f [x T A(t)x+2x T N(t)u+u T B(t)u]dt (1.1) 

t o 

Subject To: 3fc=F(t)x+G(t)u, x(t 0 )=x. 0 ' (1.2) 

Mx f » $ 0.3) 

where x is an n-vector, u is an m-vector, and $ is a 
kin -vector. The initial and final times t^ and t f are 
assumed to be known , and the matrix B(t) is assumed to be 
symmetric and positive definite on {t Q ,t f ] . It is a necessary 
condition that B(t) be at least positive semi-definite ; the case 
when B(t) is not positive definite is called the singular case 

*Numbers indicate references listed in Section 5* 
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and its solution is presented in Ref. 7. If one is only 
interested in using the computer program, then one can proceed 
immediately to the self-contained Appendix A. 

In Section 2, the solution for the problem defined by 
Eqs, (1.1)- (1.3) will be solved without recourse to optimal 
control theory methods (i.e., without using the calculus of 
variations, or the Pontryagin maximum principle). To demonstrate 
the essential features of the justification, the quantities ■ 3? 
and N(t) will be assumed to be. zero in Section 2. The case v/ith 
S? and N(t) included will be justified in Section 3 and optimal 
control theory will be utilized. It should be noted that the 
usual application of Eqs. (1.1)-(1.3) is with J and N(t) equal 
to zero, so Section 2 should be sufficient justification for 
most of the cases which arise in applications. 



2. METHOD JUSTIFICATION: WITHOUT 

OPTIMAL CONTROL THEORY 

In this section the optimal feedback control for the 
problem defined by Eqs. (1.1)-(1.3) with M=0 and N(t)=0 
will be determined, i.e,, 

Minimize: J=lx|s f x f +1 J" ; f [x^A( t)x+u T B(t)u]dt (2.1) 

t o 

Subject To: x=F( t)x+G( t)u, x(tQ)=x Q | (2.2) 

That is, we wish to determine the control 

u=f(t,x) (2.3) 

i.e., a feedback control function, which causes J (the perfor- 
mance index) to be minimized. In Section A. 2 of Appendix A 
the typical origin of the terms in Eqs. (2.1) and (2.2) is 
discussed. Also, as noted in Section 1, B(t) is assumed to 
be symmetric and positive definite (the nonsingular case) . 

We could use the calculus of variations or the maximum 
principle to solve this problem, however knowledge of varia- 
tional theory is required. Instead, we shall employ a T, trick u 
which is employed in many types of optimization analysis. 

This trick involves the introduction of an arbitrary function 
with specified continuity and differentiability properties, 
which will be chosen later to help us out. (Such a trick is 
also used when one introduces Lagrange multipliers into an 
optimization problem, i.e., they are first treated as arbitrary 



k- 


functions and then particular functional forms are chosen to 
aid in the solution of the problem.) 

Let S(t). he an arbitrary dif ferentiab3.e , symmetric 
nxn matrix function. ^ 

Property: i f ^~-^[x^S(t)x Jdt-i-[x^S( t)x] ^ S 0. (2.4) 

' J t " o 

o 

T 

Proof : The integral is an exact differential in x Sx, so: 

iJ^ f d[x T Sx3-i[x T Sx]^ f - i[x T Sx]^ f - i[x T Sx] t f ~ 0. 
x o o o o 

If Eq. (2.4) is added to Eq. (2.1), the problem is not 
changed because Eq. (2.4) is identically zero. Thus, performing 
this addition: 

J=ix^S f x f - ftx T S(t)x] t f + i i t f [x T Ax+u T Bu+ ^(x T Sx)]dt. 

Note that S f is a given nxn matrix, whereas S(t f ) is (as of now) 
just an arbitrary nxn matrix function. Combining terms outside 
the integral and differentiating under the integral results in: 

^ ( S f -S ( t f ) ) x f ■ +ix^S ( t^) Xq 

+ i T j. ‘ f [x T Ax+u‘ L Bu+x T Sx+x^Sx+x T Sx]]dt 

But, x=Ex+Gu (note this is the point where the constraints get 
into the problem) , so : 

J=^[S f -S(t f )3x f+ ix^S(t 0 )x 0+ iJ'^[x T Ax + u T Bu + xVsx + uVsx 
+x T Sx+x T SFx+x T SGu 3 dt 


or, 
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J=+2(S -S(t ) ]x +^2s(t )x +Jf J ^ f [x T (A+F T S+SH-SF)x+x' r SGu 

+u T G T Sx+u T Bu]dt (2.5) 

Note that the integrand is a quadratic form in x and u, and 
that (as of now) S is an arbitrary matrix. We shall now choose 
S(t) in such a way that the optimal control is obvious. 

If we can write the integral in J as: 


J t Kx+Lu ) ^ (Kx+Lu ) dt 


( 2 . 6 ) 


with L ^ existing, then u^-l/^Kx is the optimal control because 
the integrand is the square of (Kx+Lu), which implies zero is 
the smallest value the integrand can take and u=-L Kx causes 
the integrand to equal zero. 

Let us now expand Eq. (2.6) and equate it to the integrand 
of Eq. (2.5); this will then imply how we should choose S(t) 
to get the obvious control solution form of Eq. (2.6): 


Pt.p m P L '-P r P r P m m rn m rn. m 

J (Kx+Lu) 1 (Kx+Lu) dt- J ^ ~ [ x 1 K 1 Kx+x A Lu+u A L Kx+u 1 L 1 Lu ] dt (2.7) 

0 o 

Equating terms with the integrand of Eq. (2.5) implies 


k t k=a+f t s+s+sf 


k t l=.sg 


T T 
L i K=G 1 S 


L T L=B 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


Since B is symmetric and invertible, L is also symmetric and 
invertible, i.e., 

LL=B L 2 =B =$■ L = bL 
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By Eqs . (2.9) and (2.10) 

rp m _i T -.4. T 1 ' . . 

m'Ss^U ‘G^B ^G x 3. (2.12) 

Finally, by Eqs. (2.8) and (2.12): 

1 rp rp _JL rp rn p 

(B ^G 1 S) i (B 2 G A S)=A*F 1 S+3+SF 

or 

s+sf+f t s+a=sgb - ^b~^g t s , 

—4- 

where the symmetry of and S has been used. Then, 

s=-sf-f t s-a+sgb -1 g t s (2.13) 

Since x^ depends implicit3„y upon u, we can remove the implicit 
u-term from J by choosing: 

S(t f )*S f (2.14) 

To summarize, then, if one defines the Riccati equation : 

. s=-sf~f t s-a+sgb“ 1 g t s. 

with boundary condition: 

S(t f )*S f , 

the quantity J may be written as: 

■J=-ix^S(t„)x + J . f (B"^G T Bx+B‘"^u) T (B"^G I 'sx+B“^u)dt . (2.15) 

O O O v 


Then, the term outside of the integral is independent of u 
(since x Q is specified and S(t Q ) is well-defined by the solution 
of Eqs. (2.13) (2.14)) and the smallest possible value of the 

integral is zero, and the integral is zero if: 

B~^G T Sx+B~*u=0 


(2.16) 


• * 





Equation (2.16) defines the optimal feedback control, and 

Jsk S(t )x is the value of the performance index due to the 
0 o o 

optimal control. 

Example : Fin.: J=-]- J^(ax^+bu 2 )dt (2.17) 

Sub. to: x=u, x(o)=x o> T specified (2.18) 

Note that the system is a linear, time -invariant system, and if 
classical linear control were employed here, a linear feedback 
control with constant gains would be the typical result. The 
usual result with LQP theory is a linear feedback control with 
time-varying gains. However, by letting t f =T be large, the gains 
are approximately constant. (Constant gains may be obtained by 
choosing T= oo ; see Ref. 5*) 

Let us now compute the solution of the problem defined by 
Eqs. (2.17) and (2.18). As noted above, if T is finite, then 
time varying gains are obtained. However, we should expect the 
gain to approach a constant as T becomes large. The solution 
to this problem is defined by Eqs. (2.13)> (2.1A), and (2.16): 

u=-B -1 G T Sx (2.19) 

-1 1 

with: B = — , G=1 , and S is the solution of the Riccati 

equation: 


The solution of the Riccati equation is: 


s(t) = 


1 _ e 2'£7b’ (t-T) 
^ 3B ^ +e 2^V Ct'-T) ' ] 


( 2 . 20 ) 



Note that for T»t, S(t) m Jlib 1 . Plots of various S(t) as 
T varies are shown to the rig! 

<S(t) is basically constant 
except for a transient near 
t^sT. This behavior is also 
typical of more complicated 
time-invariant systems. 

Thus, for this problem, S(t)?^ -vfaF if T»t and the 
approximate optimal linear feedback control is (from Eq_. ( 2 . 19 )) 

u=-(*^) (1 ) > [ab* x=-|a7b > x . ( (2.21) 

This result agrees with intuition in that: (i) u is negative, 

if x is positive, which implies that the control attempts to 
drive the state x to zero; (ii) u is proportional to a/b . 

The latter result implies that if a»b, then . there is more 
weighting on the state in the performance index and the result 
is a large control value to maintain a small value of x. On 
the other hand, if a<<b, then there is more v/eighting on the 
control in the performance index and the result is a small 
control value (relative to the value of x). 






3. METHOD JUSTIFICATION : WITH 

OPTIMAL CONTROL THEORY 


In this section the optimal control problem defined by 
Eqs . (1.1) -(1.3) will be solved in general. The solution 
technique is well known, and Is similar to the developments in 
References 1 and 6. The underlying optimisation theory is 
discussed in Refs. 1, 3> 6. 

To develop the desired solution it is convenient to adjoin 
the terminal conditions, Eq. (1.3), to the performance Index, 

Eq. (1.1), with the constant Lagrange multiplier vector q, and 
the differential equations, Eq. (1.2), with the time-varying 
(in general) Lagrange multiplier vector p. Then, the augmented, 
performance index is: 


J=ix^S f x f ,+q T (Mx ;f .- 2 )+i 


f [x T Ax+2x T Nu+u T Bu 


+p T (Fx+Gu-x) ] dt . 

The Hamiltonian for this problem is 
Hfs-Jr ( x^\Ax+2x%u+u T Bu ) + p^ ( Fx+ Gu ) , 
and the resultant necessary conditions of optimality are: 


(3.1) 

(3.2) 



(3.3) 


H =0 
u 

p f =S f x f +M T q , 


(3.4) 

(3-5) 


where IL^ is considered as an (nxl) vector and is an (mxl ) 
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vector. Since the problem is nonsingular (i.e., B(t) is 
positive definite), the control, u, may be eliminated from 
the problem by Eq. (3.4), i.e., 

Vfl'WoVo (3.6) 

which implies 

u=-B~* (N T x+G T p) . (3.7) 

Then, upon substitution into Eqs. (1.2) and. (3*3), we have 

x^Fx-GB" 1 (N T x+G T p) 

p-r-Ax+NB " 1 (N T x+G T p)-F T p 

or, 

i=(F-GB" 1 N T )x-GB“ 1 G T p (3.8) 

p=(NB“ 1 N T -A)x+(NB~ 1 G T -F T )p ( 3 . 9 ) 

Since Eqs. (3.8), (3-9) > (1.3)? and (3.3) are linear in x, p, S , 

and q, it can be shown that there must exist linear relationships 
among the variables, and we introduce the unknown (for now) 
matrices Q, R, S, and V involved in these relationships 

p(t)«S(t)x(t)+R(t)q ( 3 . 10 ) 

4? =V(t)x(t)+Q(t)q. (3.11) 

(It can be shown that the resultant S(t) is symmetric , 1 and 
we shall assume this now to ease the notation.) Thus, among 
the 2 n+ 2 k -variables x, p, q, and jp there exist n+k independent 
variables, which by Eqs. (3.10) and (3.11) have been chosen to 
be x and q. Upon substitution of these relations into Eqs. (1.3) 
and (3.3) (i.e., the terminal boundary and transversal! ty 
conditions), we obtain 



or 


Mx f r=V(t f )x f +Q(t f )q 
S(t f )x f +H(t f )q=S f x f +M T q 

(K-V(t f ))x f +Q(t f )qsO (3.12) 

( s (t f )“S f )x f +(R(t f )-M T )q=:0 , ( 3 . 13 ) 

which are identities in x f and q. This implies that the 
coefficients of Eqs. (3.12) and (3-13) must vanish, and thus 

V(t f )=M , Q(t f )=0 , S(t f )=S f , R(t f )=M T . (3.14) 

Equations (3.14) define boundary conditions for the unknown 
matrices. If differential equations could be developed for the 
matrices, then the matrices could be computed by integrating the 
resultant equations backward with the boundary conditions of 
Eq. (3.14). We shall now determine such a set of differential 
equations by differentiating Eqs.' (3.10) and (3-11), substituting 
the results into Eqs. (3.8) and (3 - 9 ) , and then interpreting the 
resultant forms. 

Before we make these computations, it is instructive to 
answer the question of why Eqs. (3.10) and (3.11) were introduced 
in the first place. Our goal is a feedback control, say 

u=g(t,x,$ ) . (3.13) 

We are guaranteed the existence of relationships of the type 
assumed in Eqs. (3.10) and (3*11) by properties of linear 
differential equations (i.e., Eqs. .(3.3), (3.9)) with linear 
boundary conditions (i.e., Eqs. (1.3), (3-3)). Such relations 
are desirable because if S(t), R(t), V(t), and Q(t) can be 
determined, and if Q(t) is invertable, then the optimal feedback 



control can be computed from Eq. (3.7), l.e., 

u=-E -1 [K T x+G T Sx+G T I?( 2" 1 $ -Q _1 Vx).3 . '(3.16) 

l/o shall come back to this equation after determining the 
defining differential equations for S, R, V, and Q. 

First differentiate Eq. (3*10) and substitute the result 


into Eq. (3.9)3 i . e . , 

Sx+Sx+Rq=(MB~ 1 N T -A)x+(NB'’ 1 G T -F T )(Sx+Rq) . (3.17) 
Then upon substitution for x (from Eq. (3.8)) and upon rearrange- 
ment, we obtain 


[ S+SF+F T S+ A- ( SG+N ) B“ ' 1 ( SG+N) T ] x 

+ [R+(F T -(SG+N)B~kbE0q=O (3.18) 

Since x and q are the independent variables, Eq. (3.18) is an 
identity in x and q, and thus, the coefficients must vanish, 
which implies: 

■S=-SF-F T S-A+(SG+N)B -1 (SG+I0 T , S(t f )=S f (3.19) 

R=[-F T + (SG + lOB- 1 G T la , H(t f )«M T , (3.20) 

where the boundary conditions at t f are obtained from Eq. (3.14) . 

The equations for Q and V are obtained by differentiating 
Eq. (3.11) (noting that $ is a constant) and substituting for 
x, which gives 

r v+ v ( f-gb -1 ( n t *g t s) ) 3 X+ r q-vgb” 1 g t r] q=0 (3.21) 

This is also an identity in x and q which implies that the 
coefficients must vanish. Since the equation for V is the 
transpose of H,V(t^) r=M«R( t^) J ‘ , and S is symmetric, it follows that 
V(t),R(t) T , (3.22) 

so the. variable V(t) is eliminated. The equation for Q, with 



boundary condition from E q. (3*14) > is 

|=H T G3" 1 G T t; , Q(t f )=0 . (3.23) 

The resultant optimal control is then: 


Ur=-3 _1 rK T +G T (S-R3 1 i 1 ')]x-B"' 1 G T Rr 1 i , (3.2/*) 

where .3, R, and Q can be determined by Eqs. (3-19)? (3«20), and 

(3.23). 

Rote that Q occurs in the optimal control, Eq. (3.24), only 
in the product RQ*" 1 . This motivates one to develop differential 
equations for (S-RQ" 1 ^) and RQ” 1 (as opposed to 3 , R, and Q) , 
and it can be shown 1 that the resultant differential equations 


are exactly like the S and R equations. However, one cannot use 
these until t f - £ ( t > 0) because Q" 1 does not exist at t f 
(since Q(tJr-lO). Thus, in LQP, the 3, R, and Q equations are 
integrated backward for a small time increment, and then a 
switch-over to the direct computation ox 3-RQ R and RQ is 

made. This, of course, saves computer time. 

_i t 

Finally, it should be noted that S-RQ R can become 
unbounded. This means that the proposed problem does not possess 
an optimal solution (or a unique optimal in very special cases), 
and the time at which S-RQ R becomes unbounded is called a 
con.1u.gate point . The program, LQP, prints out the occurrence 
of a conjugate point and stops the computation. This is another 
reason for choosing S f > 0, A(t) 0, N(t)sO, and no terminal 
conditions because then one is guaranteed that a unique optimal 
control exists^ and no conjugate point can occur. 



u . 


SUW-iARY AND CONCLUSIONS 


This report contains developments ox the linear quadratic 
optimal control problem, one of which does not involve optimal 
control theory. The theory is applicable to the development of 
neighboring optimal feedback guidance gains, and is useful as a 
tool for synthesizing feedback control laws in general. A 
computer program- which requires only the pertinent matrices of 
the linear quadratic problem is described in Appendix A, which 
also serves as a self-contained User’s Guide. 

Knowledge oi optimal control theory is not necessary to use 
the computer program or to understand the development of the 
expression for the optimal feedback control (see Section 2). 
Thus, Section 2 and Appendix A may be learned in a relatively 

snort period of time without any background in optimization 
theory. 

The relationships between classical feedback control design 
and linear quadratic optimal control design were presented in 
a number of lectures to KASA-JSC and contractor personnel in 
July- August 1974 by 17. F. Powers. Lecture Notes were handed out 
at the lectures and are available upon request, from Modern 
Systems Analysis, Inc. 


1/ 4 
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APPENDIX A 


USSR’S GUIDE FOB 


iiQP is a subroutine which solves the following optimization 
problem, which does not require iteration: 


Minimize : 


x T A(t)x+2x x K (t)u-Ki T 3(t)u]dt 


Tf f 


(A . 1 ) 


o 


Subject To: x«F( t)x+G( t)u , x(t )=sx (A. 2) 

vJ w 

Ksc f « f (A.3) 

where x-n-vector, u-m-vector, $ =p~vector, and t Q and t f are 
specified. The notation of Fqs, (A,1)~(A.3) is that of Ref*. 1 . 

The user need only supply a ’'MAIN 1 ' subroutine which defines the 
parameters of the problem and calls LQP. If any of the matrices 
A, IT, B, F, or G are time-varying, then a second subroutine 
TIMVAR, which defines the time-varying matrices, must be supplied, 
also. Since LQP employs the numerical integration scheme DVDQ 
(which is a variable-stepsize , variable -order integrator; see 
Ref, 2), it is recommended that the time -varying matrices in 
TIKVAF be represented by cubic splines. 

A.l Basic Floy/ Of The Algorithm 


As shown in Section 3, the solution of the optimal control 
problem defined above is: 


1 rpm -tm i rn 1 . 

P * r TiT - . f - ( C T)A“ 'ti 1 ^ 1 -n~ > n 1 Vp 

U=T-iJ J-! “VI £•£>£ n ) jX— £J Cr aAv X 


(A. A) 


where 
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S^-SF-^S-A+CSG+I'OB -1 (SG+N) T , S(t f )=S 
R=[-V T +(SG-tJ?'3B" 1 6 T ]i , R(t M# 


« fj 1 1 f T 1 

P r*- L r*T> • n — ■> 


Q=F GB G 1 . , 0 ( t ,r> ) “ 0 


(A. 6) 
(A. ?) 


If, as in most applications, Eq.. (A. 3) is not present, then the 


solution is defined by: 

_1 rp m 

u*-B l rN A +G*S]x , 


(A. 8) 


where S is still defined by Eq. (A. 3). Since the latter problem 
requires much less integration and logical operations, it is 
advantageous to model the control problem without an Eq. (A *3) 

(if possible) and a flag exists in the program for this purpose 
(I FLAG! ) . 

The flow of the computations is as shown in Figure A.1, 
i.e., the values for and Q(t f ) are defined, numerical 
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integration proceeds backward to t Q , and then the optimal state 
and control are defined by a forward integration. 

4*2 Select ion Of Weighting Matrices 

In the optimization problem of Eqs. (A,1)-(A.3), the 
matrices E(t), G(t), M, and 4* are defined by the process, 
e.g., F and G typically result from linearization about a nominal 
trajectory . If pure neighboring optimal guidance is to be used, 
then the matrices in Eq. (1) are also well-defined (e.g., see 
Refs. 1, 3)* However, most applications will probably require 
the specification of the weighting matrices S f3 A(t), H(t), and 
B(t) by the guidance or control designer. In this section a 
,r rule- of- thumb’' for weighting matrix selection which has proved 
useful in a number of applications (Ref. 1, 4) will be presented. 

To get started on a design, assume N(t)=0, i.e., no mixed 
state-control terms in Sq. (A.1). In most cases one will not 
have to employ a nonzero H-matrix at any Unite in the design. The 
only remaining matrices are S f , A, and B which weight terminal 
state values, state trajectory values, and control values, 
respectively , If Eq. (A. 2) results from linearization about a 
nominal trajectory (the usual case) , then x and u actually, 
represent deviations from the nominal. In such a case, one usually 
has some idea of the tolerable deviations for each variable. Thus, 
assume it is desired that: 


IV- 

|x. j (1=1,. . . ,n) 
1 1 f l 

(A. 9) 

K (t) t 

£ |* ± (t)| 

(1=1 , . . . ,n) 

(A. 10) 

| u i (t) l 

< (Vt)| 

U = 1 > * * * jKl) 

( A . 1 1 ) 
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i.e., the maximum deviation from the nominal value associated 
with x.j is _+ and so on. Then, the smaller the value, of 

TT (e.g., x^) , the larger the weighting of () (e.g., „) 

should be in Eq. (A.T), and vice versa. A choice which satisfies 
this criterion is: 



B( t) ss 




(A. 14) 
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Typically the x\ (t) and IL (t) values are constants; however in 
Shuttle reentry one may wish. to change the weighting matrices 
from one flight phase to another (e.g., from the constant dra.g 
phase to the equilibrium glide phase and so on) . One can then 
compute the resultant optimal control with LQP, and check to 


see if the resultant feedback control meets all specifications. 

If not, the weighting matrices should be modified, and, of course, 
the modifications are problem dependent. In any case, Eqs. (A.12)- 
(A.Ul) give, at least, a well-defined start to the feedback gain 
design process. 

Finally, to save computer time, B”* 1 wall be supplied to 
the program instead of B. Usually is easily calculated 
beforehand (if not, the computation can be done in TIMVAR) . 


A. 3 LQP Argument List 

In this section the variables employed in LQP will be listed 
along with their type (integer or double precision), dimension, 
and identification with the variables in Eqs. (A.1)-(A.3). A 
T: MAIN n subroutine (to be discussed in the next section) is to be 
supplied by the user, and a call to LQP is made from M&IN. The 
CALL-statement is: 

CALL LQP ( N , K, IP , I FLAG 1 , I FLAG2 , TI , TF , EP ,SF,A, DK , BINV , F . G . 

DM, PCI , X , SRQX, DSRQX, S,R,Q,U, XDOT, DT, KQ, YN , DUM1 , 

DUM2 , DUM3 , DUMA , DUM3 , DUM6 , DUM? , DUMB , DUM9 , DUM1 0 , DUM1 1 , 

XL , XLDQ T , DUM1 2 ) ' 


The variables in this call are defined as follows , where I^integer, 
Dndouble precision, 3* single precision, and 


K=n(n+1)/2 + np*p(p*1 )/2+n 



Program 

Variable 


Problem 

Ya.ri3.ble 


m 


I FLAG 1 

IFLAG2 

TI 

TF 

SP 

SF 

A 

DN 

BINV 

F 

G 

DM 

PSI 

X . 


O.L\'C'A 


DSPQX 


o 


S. 

A 


M 


4 

x 


-1 


Q 


u 


v'TV'-Ui 


n 

V.' 

U 

w 

X 


DT 


Variable 


I 

T 


J. 


I 

D 

D 

S 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 

D 


■ Dimension 
scalar 
scalar 
scalar 
scalar 
scalar 
scalar 
scalar 
scalar 
N*(N+1 )/2 
'N*(N+1)/2 
N*M 

M * ( M + 1 )/2 

N*N ' . 

N*M 

N*IP(If p=sO,.DIM=l ) 
IP(If p-0,DIFM1) 

N 

K 

K 

N*(N+ 1)/2 

N*IP(If p=0,DIM=:1) 

IF*(IP+1)/2(If p-0, 
DIIi= 1 ) 

if 

1 '! 

07,10 


REPRODUCIBILITY 
ORIGINAL PAGE IS 


OF TH E 
POOR 



Program 

Variable 

Problem 

Variable 

Variable 

Type 

Dimension 

K Q 

— 

T 

j. 

- ‘ 

fl'-i 

— 

D 

TV 

DU Ml 


D 

!I*N 

DUM2 

-- 

D 

IT*K 

DUI43 

-- 

D 


DUMZf 

— 

D 

N*M 

.DUI45 

— 

D 

N*M • 

DUM6 

— 

D 

p^O , DIM- 1 ) 

DUH7 

— 

D 

L*I P (If >psO , DIM- 1 ) 

DUM8 

— 

D 

N*IP(If psO, DIMsl ) 

DUM9 

— 

D 

N*IP(lf p— 0 j Dlilrr 1 ) 

DU Ml 0 

— 

D 

M 

DUM1 1 

— 

D 

M*(M+l)/2 

XL 

p(see Eq.3.9) 

D 

N 

XLDOT 

p(see Eq. 3 . 9 ) 

. D 

II • 

DUM12 

— 

D 

N 

The 

variables above which 

are not problem 

variables are 

described 

below. Except for IP 

'LAG! , IFLAG2, and EP, these 


variables are LQP and DVDQ "working variable s K which are of no 
concern, to the user except for DIMENSION statements (and the 
dimensions are well-defined in the list above) . 

I FLAG 1 : flag set by user indicating presence or absence of 

t e rmin al condi tion s ; -0 if terminal c ondi tions pre sen t , 
and r; 1 if no terminal conditions. 



IFLAG2; flag set by user indicating time- variability of' 

matrices; -0 if A, 3, DTI, F, G are time -invariant, 
and rrl if at least one of these matrices is time- 
varying (in which case the user must supply a subroutine 
TIKVAR) . 

BP: absolute local error indicator for the numerical integration 

scheme DVDQ; this parameter is problem dependent, but a 
safe initial choice is 1 .2-5* (Gee Ref. 2 for a more 
thorough description of EP.) 

SBQX: contains the vector being integrated by DVDQ. If terminal 

conditions present , SRQX contains S, R, Q (expressed in 
vector form) going backward and 5, II, Q, x going forward. 
If no terminal conditions, SRQX contains £ going backward 
and x, p going forward. 

DSROX: contains the time derivative of SRQX. 

DT: storage required for DVDQ. 

KQ : storage required for DVDQ, . 

YN: storage required for DVDQ . ; 

DUH1 through DUI-112: dummy storage required for LQP matrix 

manipulations. 

A . A Matrix To Vector Conversions 

Sven though Eqs. (A.I)-(Ao) are written in matrix form, the 
computer program operates in a vector mode. (The only matrix 

o 

dimension is for DT, which is part of the integrator, DVDQ). 

Thus, all matrices must be converted to vectors, and since some 
of the matrices are symmetric, considerable savings can be gained 
b 3 r distinguishing between general and symmetric matrices. 



number of subroutines 


A number oi subroutines from Kef. 8 are employed in LQP 
■to perform the various matrix manipulations. These subroutines 
assume that the matrices have been converted to vectors column- 
by~column . That is, consider the three-by-three matrix A: 


a 11 

a 12 

a 13 

a 21 

a 22 

a 23 

a 31 

a 32 

a 53 

matrix, 

then 

it wi 


(A. 15) 


co lunm-by- column , i . 0 . } 

T 

A(9) - r a-| -j ^21 ^3 ^ *j o ^22 ^"^2 a i3 ^23 a 33 ^ * ( A » 1 o ) 

(General Matrix Format) 

If A is a symmetric matrix, then it will be converted to a 
6- vector column-by-column of the upper triangular portion of the 
matrix, i.e., 


A(6Ma n a 12 
(Symmetric 


a 22 a l3 a 23 a 33 ] 
Matrix Format) 


(A. 17) 


The various matrices are printed out in the same manner (i.e., as 
n or n(n-M )/2 vectors in the format of Sqs. (A. 16) or (A. 17) > 
respectively) . 

A. 5 Example Problems 

A number ox simple examples will be presented in this section 
to illustrate the setup of MAIN and TIMVAP, . typical printout, 
and the output for a problem path a conjugate point. 



Sxam-ole A.1: Let x be a 3-vector and u a 2-vector 


Itin.im.se : 




(2x?+u^+Up) dt 


( n 1 ON 
V. -ri. . 1 Kj J 


Subject To: 


x* 


0 1 0 
OOO 
0 0 1 


0 0 

1 0 
0 1 


u 


(A. 19) 


x 1 (0) ~1 

, x 2 (0) 

=2 , 

1 0 0 


0 

0 1 0 

X(1) r 

0 

0 0 1 


1 

- _ 


_ tm 


3U(0)=0 



(A. 20) 


(A, 21) 


Equation (A. 1 8) corresponds to ICq. (A.1); Eqs. (A. 19) and (A. 20) 
to Sq. (A. 2); and Eq. (A. 21) to Eq. (A. 3) • Since terminal 
conditions are present, IFLAG1=Q; and since all of the matrices 
are time-invariant, IFLAG2-Q. A value of EP^I .0x1 will be 
used for the absolute local error control in the integrator. A 
typical MAILT, subroutine ( the only information . required by the 
user) is shown in Figure A. 2. The development of MAIN involves 
the development of a well-defined DIMENSION- statement, data 
input, and a well-defined CALL to LQP. (Also, note the HEAL EP 
statement because EE must be single precision to avoid difficultie 


on UNI VAC computers.) 

The program begins the backward integration of S, 

at tsl , and the printout is shorn in Figure A. 3. Since 
<- T T ~ - 1 

Sr*.S-KQ R and are needed to define the optimal 


Q, and- R. 
only 

feedback 
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gains, a switch to the .3 and V/ system is made at t-0,93* 

(Such a switch is possible only if Q" 1 (0-95) exists; if Q _1 (0.95) 
does not exist, then the program will stop, usually indicating 

/V» 

an abnormal problem. ) From t=?0.93 to t~.Q, the 3 and VJ matrices 
are printed out in the format of Sqs. (A. 17) and (A. 16), 
respectively, since S is symmetric and A' is a general matrix. 

The first few forward integrations from t».0 are shovrn in 
Figure A . 4 , and. the values of the resultant optimal state and 
control are added to the printout. One can compare the backward 
and forward values of S and 17 to aid in the choice of a value 
for FP which gives the desired accuracy. 

Example A. 2: Let x and u be scalars. 


Minimize : 


- J 


23L 

2 


0 


u 2 dt 


Subject To: xsx+u .■ 

x(0) =0 , x ( 3 tt/2 ) ~ 1 . 0 

This optimal control problem possesses a conjugate point at 
t= TT/2, which implies that there does not exist an optimal 
feedback control on the interval [0, 3 tt/2 ] . The program detects 
the possibility of a conjugate point when the numerical integra- 
tion scheme begins to decrease the stepsi&e to a very small value 

IV 

(which' is necessary to get accurate values of 3 and W since 
(si — > <>d ) . The output of the program denoting that this 

behavior is occurring at t~ tt/ 2 is shown in Figure A. 3. 





Let x and u be scalars. 


Minimise : 



Sub j ec t To : x- t L "x* tu , x ( 0 ) - 1 

This optimal control problem possesses time variable dynamics , 
and, thus, the subroutine TIMVAR is required. A typical T1MVAR 
subroutine for this problem, is shown in Figure A. 6. 
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IMPLICIT ft '£ A L * 8 ( A ~ H » 0 - t ) 

OiAnsION SF'(i»',*(6t,0N«6),8lNV(3»,F«9l ! • ® «'* * * “7 1 S t I W t 2 H-i-I K »■€ 
v>, i 7U J nU H ! < 9 ) , 0U«2 (<?) ', DUft3 t 9 > t DUNH J 6 > » OUiib * 6 1 * ol,1 ° ' ° 1 

0 U H a (-9-! ■% U M 9 ( 9 ) o u H 1 0 ( 2 J ■. ou H 1 1 ( 3 I , VUS-'i , X L “ rW J DWM CV> 


N - 3 ■•" 

rt = 2 

•I P = 3 - - • 

If LAG l !a O 
! F L A G 2 = 0 
T l “0 a DQ 

TF»1 •DO- 
£P=*i •> E -5 
0 0- i ■!»!'• 6--; 
SFU) =0,00 

0 0 2 ■ | * l i • & ' 

A 1 I ) ■« G - 0 0 


■ A { 6 } *2*00 

0 o 3 1 - 1 * 6 
ON 5 i ) »b* D 0 
3 l h V ( 1 > * 1 * 0 0 
a i N V { 2 1 * o ♦ o c 
3 l n V { 3 ) * l a 0 0 
DO H 1-1*9 - 
f { i ) -0*0Q 

F t *t ) ” 1 «D.Q 

f ? 7 ) - 1 »oo 
005 1-1*6 " 

GU)®0«00 
- G i 2)” l *00" 

G 1 6 ) » 1 * 00 
DO 6 1-1 '» 7 " 
DM ( 1 ) -0.00 
°mu ) = r. D o 
OfH 5 ) * 1 • DO 
••• - DM C9")»T* DO 
p's l ( l ) 3 0 • 00 

P5 I ( 2 ) * C • 0 0 
PS I ( 3 ) a 1 * DQ 

x m a i t do * 

X i 2 » =2*00 




* 



LQp?N*M*lP*IFl> AGl,lFuAG2$ 
X * 5 R <3 X I 0 S R Q X . s * ft » 0 * U t X D 0 T 

- n 1 1 u o imimq .nUMlfliOUM 


H ,Tf tEPiSF.AiUMiBlNV.F 

. - _ is rt v m n»Mi OUM 


END 


Figure A. 2 
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->l MAIN Subroutine Eor Example A.1 



1 ® coo 




* ^ 


backward integration begins at r~ 


T “ roOQQ 


s mat r t x : 

• oooooo 

0 OOOOOO 

*00000 0- 

R matrix: 

l » cocooo 

® OOOOOO 

. QCOOQO 

R matrix: 

•OQOOCQ 

l.OQOCOO 


u matrix: 

•OQQOOQ 

.oocooo 

« OQOOOO 


T= « 95 0" - - - - - - - •; — - - •' _ 

s matrix: 

.000000 

. OOOQQO 

•OOOOOO 

H matrix; 

1,000000 

.050000 

.OOOOOO 

R matrix: 

.OoOoOO 

1.048559 


q matrix: 

-.000042 

'*‘'001250 

“.05G000 


SUCCESSFUL SWITCH 

TO T ILOE SYSTEM aT 

T® .9 50 

T * *9 50 








— s M jlde) : 

-9641 7.826274 

-24 10.445657 

80.26 1 1-4 1 

----- W I-T-I lD£ )■'•••• 

- -*9 64 l 7.826-27-4— 

“24 l 0.-445657 

.00 DO 00- 


• uuyuuu 

- 1 y » v / bu^2 ' 


T » *902 

Sl'TJLOEJ l 

1 2 9 54 b 5 l 8*704- - 

6 3 1.53279 7 - ; - 

41.043634 • 

• -wf tilde) ; 

-1 295 4,5i8?04 - 

-6 31 .5 32 7 97- — 

- — .OOOOOO- 

ii I T ♦ L U £ I • 

• OCOG'OO 

- i C * 2078 22 



Figure A. 3 Backward Integration Printout For 

Example A.1, With Switch At t=0.95. • 


<- o 
*> <■ 

• e r 


#< 
«■ r 

* «• 




. #o*— 
t** 




0 ®G^ C ®II JTr 




FORWARD ~INT£G«AT i OH STASIS 


T- »QQD 

-» J * * W ' ■ " "" 

SITILD6): - -12*000006 • 6oOGQOQ3 : - 

W.U-U.0E4 *-•■ — ....^1.2,000006 -6,000003- - 

W-CT-l-tO-E ll~ : *000000- --*&3267-7 - 

~ — state:- ' — - — i • oooooo---- — • 2 • cooooo ; — - 

— CON-TROLl •» 1 H *000002 *6 32 6-7-7 


T 3 . * 0*1 7 

- - 


s niuoe) 5 - — 

— 13 *868301 - 

- 6*613351- 

Vii Tl^og-H - 

-13 * 8B63Q1 

-6*613351 

h l'T i LOE l : 

bOOOOQQ 

• -* 6 90923 ------ 

STATES- - 

1*079635 

1 0 3 6 2 0 7 5 -• - 

control: 

** 12 * 8 6 OoO 6 

• *604733 

Figure 

A , h- Forward 

For Exar 

Integration Printout 
nrle A.1 . 


H. 000Q02 - 
- *OOOGGCL 

*0000 00 *: 


4*199477 
• OOODOO 

0 0 3 0 0 8 6 



I-: -32 

- J , 4 4 •> v V 3 

reproducibility of thf 
original page is poor 1 

S (TILDE) : 

-15.361161 ' 


w (TILDE) : 

-16.391674 



i r i » t * » !!!!!!!!!!!!!!!!!!! 

- * 

! DVDQ SUSPENDS EXECUTION WITH IFLAG = 7 , i 

! ! 

i If IF LA 3 = 7 (MINiaUM STEPS IZE EXCEEDED)!, * 
I * 

J LI K ELY A CONJUGATE POINT EXISTS AT ! 

t * 

.! T = -1. 57342 3 

: ■ . 

!! I! !!!! I !!!!!'.!!!! H !!! J !!!!!! 1 !!! J H !!!!’ !? 

Figure A. 5 Printout Denoting The Occurrence Of A 
Conjugate Point In Example A. 2. 


SUBROUT INP TIMVAR (T , N , M, A , ON, B INV f F ,G) 
IMPLICIT REAL*8 (A-H,G-Z ) 

DIMENSION A ( I) ,ON(l) , BINV( I) , FCl ) ,CU ) 


F ( l) = T*T 
G (1 ) =' T 
RETURN 
END 


A Typical TIMVAP Subroutine' For 
Example A,;, 


Figure A. 6 



